##################################################
## Author: Stephanie L. DeMora, Christian Lindke, Sean Long, Jennifer L. Merolla, and Maricruz Ariana Osorio
## Project: The Effect of the Political Environment on White Women’s Political Ambition
## Purpose: Replication File for Publication: In-Text Tables & Figures
## Date: May 9, 2023
##################################################

# Set working directory to wherever you saved the files:
setwd("Replication Files/")

# Load in all packages
Packages <- c("here", "car", "survey", "margins", "interplot","aod","dplyr","foreign","nnet","ggplot2","reshape2","Rmisc",
              "tidyverse", "broom","dotwhisker","janitor","qualtRics","jtools","mediation","gtools","gridExtra", "huxtable", "DT", "ggsci", "hrbrthemes", "nnet", "haven", "estimatr", "ivreg", "readr", "nlme", "lubridate", "readxl", "tseries", "forecast")
lapply(Packages, library, character.only = TRUE)

# Set plot theme:
steph_theme_gray <- list(
  hrbrthemes::theme_ipsum_rc(base_size = 12,
                             grid= "",
                             plot_title_size = 20,axis_text_size = 18),
  ggsci::scale_color_material("grey"),
  ggsci::scale_fill_material("grey"),
  theme(legend.position = "top",
        plot.caption = element_text(hjust = 0, face= "italic"), #Default is hjust=1
        plot.title.position = "plot", #NEW parameter. Apply for subtitle too.
        plot.caption.position =  "plot"))

scale_norm <- function(x) (x-min(x, na.rm=T))/(max(x, na.rm=T)-min(x, na.rm=T))
margins_summary <- function(model, ..., level = 0.90, by_factor = TRUE) {
  mar <- margins(model, ...)
  summary(mar, level = 0.90, by_factor = by_factor)
}

# ---------------------------------------------------------------------------------------------------------------

# Load TESS/NORC Data 
df_TESS <- readRDS("NORC_Data.rds")

# ---------------------------------------------------------------------------------------------------------------
# ---------------------------------------------------------------------------------------------------------------

# Demographics for Appendix Table A1.

df_TESS %>%
  tabyl(AGE7, sort = TRUE)

df_TESS %>%
  tabyl(GENDER, sort = TRUE)

df_TESS %>%
  tabyl(RACETHNICITY, sort = TRUE)

df_TESS %>%
  tabyl(AGE, sort = TRUE)

df_TESS %>%
  tabyl(PARTISANSHIP, sort = TRUE) %>%
  select(-valid_percent)

df_TESS %>%
  tabyl(EDUC, sort = TRUE)

# ---------------------------------------------------------------------------------------------------------------

# Scales and Alphas on Page 19
# Fear
x <- cbind(df_TESS$M1_3, df_TESS$M1_6)
psych::alpha(as.data.frame(x))

# Enthusiasm:
x <- cbind(df_TESS$M1_1, df_TESS$M1_4)
psych::alpha(as.data.frame(x))

# Anger:
x <- cbind(df_TESS$M1_2, df_TESS$M1_5)
psych::alpha(as.data.frame(x))

# ---------------------------------------------------------------------------------------------------------------

# Anger Averages on pages 19-20:
summarySE(df_TESS, measurevar="anger", groupvars=c("TREATMENT_sort"), conf.interval = .90)
t.test(x = df_TESS$anger[df_TESS$TREATMENT_sort == "Trump"], y = df_TESS$anger[df_TESS$TREATMENT_sort == "Alexa"], alternative = "greater")
t.test(x = df_TESS$anger[df_TESS$TREATMENT_sort == "Clinton"], y = df_TESS$anger[df_TESS$TREATMENT_sort == "Alexa"])
t.test(x = df_TESS$anger[df_TESS$TREATMENT_sort == "March"], y = df_TESS$anger[df_TESS$TREATMENT_sort == "Alexa"])
t.test(x = df_TESS$anger[df_TESS$TREATMENT_sort == "MeToo"], y = df_TESS$anger[df_TESS$TREATMENT_sort == "Alexa"])

# ---------------------------------------------------------------------------------------------------------------

# Figure 1: Mean Level of Anger by Experimental Condition in NORC Sample

df_TESS1 <- subset(df_TESS,df_TESS$anger != 'NA')
tgc <- summarySE(df_TESS1, measurevar="anger", groupvars=c("TREATMENT_sort"), conf.interval = .90)
tgc$low = tgc$anger-tgc$ci
tgc$high = tgc$anger+tgc$ci
tgc$sort <- c(5,1,2,3,4)
tgc %>%
  arrange(sort) -> tgc
ggplot(data=tgc, aes(x=reorder(TREATMENT_sort, sort), y=anger)) +
  geom_bar(position = "dodge", stat="identity") +
  geom_errorbar(position = "dodge",aes(ymin=low, ymax=high)) + 
  steph_theme_gray + 
  labs(title="Mean Levels of Anger", 
       x="Treatment Condition", y = "Average level of anger", caption = "Fall 2019 NORC Survey",subtitle = "By Experimental Condition", fill = "Gender:") + 
  theme(legend.position = "none")

# ---------------------------------------------------------------------------------------------------------------
# Table A2
df_TESS$party <- car::recode(df_TESS$PARTYID7, "1 = 1; 2 = 1; 3 = 1; 5 = 0; 6 = 0; 7 = 0; 4 = NA")
Table_A2_A <- lm(anger ~ TREATMENT, data = df_TESS)
Table_A2_B <- lm(anger ~ TREATMENT*party, data = df_TESS1) # Run Linear Model examining interaction of treatment and party
export_summs(Table_A2_A,Table_A2_B, note= "+++ p < 0.001; ++ p < 0.01; + p < 0.05. one-tailed", stars = c(`+++` = 0.002, `++` = 0.02, `+` = 0.1)) #Set stars to double 0.05, 0.01, and 0.001 to show significance for one-tailed test rather than 2, note in the notes section.
export_summs(Table_A2_A,Table_A2_B,error_format = "({std.error})") # two-tailed

# ---------------------------------------------------------------------------------------------------------------

# Figure 2: Change in Anger Moving from the Control Group to Conditions by Party  
df_TESS1 <- subset(df_TESS, df_TESS$party != "NA") # Create Data Frame for Analysis
anger1 <- lm(anger ~ TREATMENT*party, data = df_TESS1) # Run Linear Model examining interaction of treatment and party
# Generate a summary of the marginal effects.
margins::margins_summary(anger1, at = list(party = range(df_TESS1$party)))[3:10,] # The [3:10,] limits reported margins to actual treaments.
# Generate a Dataframe using the marginal effects.Since this data set is Educated Women, we don't need to subdivide data frame
# by gender.
data <- data.frame("model" = c("Republican Women","Democratic Women","Republican Women","Democratic Women","Republican Women","Democratic Women", "Republican Women","Democratic Women"), 
                   "estimate" = c(-0.00226, 0.00819, 0.0467, -0.0128,0.0012,-0.0044,0.1021,0.5235), 
                   "conf.low" = c(-0.051, -0.0284, -0.0008, -0.0495,-0.0463,-0.0410,0.0516,0.4862),
                   "conf.high" = c(0.0465, 0.0448, 0.0942,0.0239,0.0487,0.0321,0.1525,0.5608), 
                   "sort" = c(1,2,3,4,5,6,7,8),
                   "term"=c("Clinton", "Clinton", "Marches","Marches", "#MeToo", "#MeToo", "Trump","Trump"))
# Plot a dot and whisker graph for the marginal effects of Anger in Black and White 
dwplot(data) +
  theme_bw() +
  ggtitle("Treatment Effect on Anger") +
  steph_theme_gray +
  theme(legend.justification=c(-1.0, .993), legend.position=c(.02, .99),
        legend.title = element_blank(), legend.background =
          element_rect(color="gray90")) + 
  xlab("") +
  scale_colour_grey(start = .3, end = .7,
                    name = "model") +
  labs(caption = "Fall 2019 NORC Survey") +
  geom_vline(xintercept = 0, colour = "grey60", linetype = 2) +
  expand_limits(x = 0.7)  # Add this line to include 0.5 in the x-axis

# ---------------------------------------------------------------------------------------------------------------

# Figure 3. Path Diagram of Mediation Analysis & Table A3
df_TESS %>%
  filter(!is.na(Ambition_norm)) -> df_TESS
# Create the dummies for consistency:
df_TESS$CLINTON <- car::recode(df_TESS$TREATMENT, "'Alexa'= 0; 'Trump'=0; 'MeToo'=0; 'March'=0; 'Clinton'= 1")
df_TESS$TRUMP <- car::recode(df_TESS$TREATMENT, "'Alexa'= 0; 'Trump'=1; 'MeToo'=0; 'March'=0; 'Clinton'= 0")
df_TESS$METOO <- car::recode(df_TESS$TREATMENT, "'Alexa'= 0; 'Trump'=0; 'MeToo'=1; 'March'=0; 'Clinton'= 0")
df_TESS$MARCH <- car::recode(df_TESS$TREATMENT, "'Alexa'= 0; 'Trump'=0; 'MeToo'=0; 'March'=1; 'Clinton'= 0")
df_TESS$Ambition <- df_TESS$Ambition_norm
mix_mod <- '
Ambition ~ b1 * fear + b2 * anger + b3 * hope + c1 * TRUMP + c2 * CLINTON + c3 * MARCH + c4 * METOO
fear ~ a1 * TRUMP + a99 * CLINTON + a98 * MARCH + a97 * METOO
anger ~ a2 * TRUMP + a299 * CLINTON + a298 * MARCH + a297 * METOO
hope ~ a3 * TRUMP + a399 * CLINTON + a398 * MARCH + a397 * METOO
indirTRUMPfear := a1 * b1
indirTRUMPanger := a2 * b2
indirTRUMPhope := a3 * b3
indirCLINTONfear := a99 * b1
indirCLINTONanger := a299 * b2
indirCLINTONhope := a399 * b3
indirMARCHfear := a98 * b1
indirMARCHanger := a298 * b2
indirMARCHhope := a398 * b3
indirMETOOfear := a97 * b1
indirMETOOanger := a297 * b2
indirMETOOhope := a397 * b3
# total effect
totalTrump := c1 + (a1 * b1) + (a2 * b2) + (a3 * b3)
totalClinton := c2 + (a1 * b1) + (a2 * b2) + (a3 * b3)
totalMarch := c3 + (a1 * b1) + (a2 * b2) + (a3 * b3)
totalMeToo := c4 + (a1 * b1) + (a2 * b2) + (a3 * b3)
'
fit <- lavaan:::sem(mix_mod, data=df_TESS, std.lv=TRUE)
summary(fit)
options(scipen=999)
par(mfrow = c(2,2))
obj <- semPlot:::semPlotModel(fit)
obj@Pars %>%
  dplyr::filter(lhs == "TRUMP" | lhs == "Ambition" | lhs == "fear" | lhs == "anger" | lhs == "hope") %>%
  dplyr::filter(rhs != "CLINTON") %>%
  dplyr::filter(rhs != "METOO") %>%
  dplyr::filter(rhs != "MARCH") -> obj@Pars
obj@Pars$rhs <- car::recode(obj@Pars$rhs, "'hope' = 'enthus'")
obj@Vars$name <- car::recode(obj@Vars$name, "'hope' = 'enthus'")
obj@Pars$lhs <- car::recode(obj@Pars$lhs, "'hope' = 'enthus'")
l = round(obj@Pars$est[1:7],digits = 2)
star = c("*","*","*","","*","*","*")
lstar <- paste(l,star,sep = "")
summary(fit)
semPlot::semPaths(obj, "est",edge.color = "#000000", fade = F, residuals = F, intercepts = F,layout = "tree",
                  sizeMan = 12.5, sizeInt = 35, sizeLat = 25, edge.label.cex=1.4, curve = 1.25,curvature = 2,
                  nCharNodes = 8, style = "ram", intStyle = "multi", title = FALSE, title.cex = 28,rotation = 1 
                  , edgeLabels = lstar
)
obj <- semPlot:::semPlotModel(fit)
obj@Pars %>%
  dplyr::filter(lhs == "CLINTON" | lhs == "Ambition" | lhs == "fear" | lhs == "anger" | lhs == "hope") %>%
  dplyr::filter(rhs != "TRUMP") %>%
  dplyr::filter(rhs != "METOO") %>%
  dplyr::filter(rhs != "MARCH") -> obj@Pars
obj@Pars$rhs <- car::recode(obj@Pars$rhs, "'hope' = 'enthus'")
obj@Vars$name <- car::recode(obj@Vars$name, "'hope' = 'enthus'")
obj@Pars$lhs <- car::recode(obj@Pars$lhs, "'hope' = 'enthus'")
l = round(obj@Pars$est[1:7],digits = 2)
star = c("*","*","*","","*","","*")
lstar <- paste(l,star,sep = "")
summary(fit)
semPlot::semPaths(obj, "est",edge.color = "#000000", fade = F, residuals = F, intercepts = F,layout = "tree",
                  sizeMan = 12.5, sizeInt = 35, sizeLat = 25, edge.label.cex=1.4, curve = 1.25,curvature = 2,
                  nCharNodes = 8, style = "ram", intStyle = "multi", title = FALSE, title.cex = 28,rotation = 1 
                  , edgeLabels = lstar
)
obj <- semPlot:::semPlotModel(fit)
obj@Pars %>%
  dplyr::filter(lhs == "MARCH" | lhs == "Ambition" | lhs == "fear" | lhs == "anger" | lhs == "hope") %>%
  dplyr::filter(rhs != "TRUMP") %>%
  dplyr::filter(rhs != "METOO") %>%
  dplyr::filter(rhs != "CLINTON") -> obj@Pars
obj@Pars$rhs <- car::recode(obj@Pars$rhs, "'hope' = 'enthus'")
obj@Vars$name <- car::recode(obj@Vars$name, "'hope' = 'enthus'")
obj@Pars$lhs <- car::recode(obj@Pars$lhs, "'hope' = 'enthus'")
l = round(obj@Pars$est[1:7],digits = 2)
star = c("*","*","*","","*","","*")
lstar <- paste(l,star,sep = "")
summary(fit)
semPlot::semPaths(obj, "est",edge.color = "#000000", fade = F, residuals = F, intercepts = F,layout = "tree",
                  sizeMan = 12.5, sizeInt = 35, sizeLat = 25, edge.label.cex=1.4, curve = 1.25,curvature = 2,
                  nCharNodes = 8, style = "ram", intStyle = "multi", title = FALSE, title.cex = 28,rotation = 1 
                  , edgeLabels = lstar
)
obj <- semPlot:::semPlotModel(fit)
obj@Pars %>%
  dplyr::filter(lhs == "METOO" | lhs == "Ambition" | lhs == "fear" | lhs == "anger" | lhs == "hope") %>%
  dplyr::filter(rhs != "TRUMP") %>%
  dplyr::filter(rhs != "MARCH") %>%
  dplyr::filter(rhs != "CLINTON") -> obj@Pars
obj@Pars$rhs <- car::recode(obj@Pars$rhs, "'hope' = 'enthus'")
obj@Vars$name <- car::recode(obj@Vars$name, "'hope' = 'enthus'")
obj@Pars$lhs <- car::recode(obj@Pars$lhs, "'hope' = 'enthus'")
l = round(obj@Pars$est[1:7],digits = 2)
star = c("*","*","*","","*","","*")
lstar <- paste(l,star,sep = "")
summary(fit)
semPlot::semPaths(obj, "est",edge.color = "#000000", fade = F, residuals = F, intercepts = F,layout = "tree",
                  sizeMan = 12.5, sizeInt = 35, sizeLat = 25, edge.label.cex=1.4, curve = 1.25,curvature = 2,
                  nCharNodes = 8, style = "ram", intStyle = "multi", title = FALSE, title.cex = 28,rotation = 1 
                  #, edgeLabels = lstar
)

# ---------------------------------------------------------------------------------------------------------------
# ---------------------------------------------------------------------------------------------------------------

# Imai mediation for Footnotes starting on Page 21:

##################### HOPE #####################

# Clinton:
df_TESS %>%
  filter(!is.na(hope), !is.na(Ambition_norm)) -> df_TESS1
model.0 <- lm(Ambition_norm ~ TREATMENT, df_TESS1)
model.M <- lm(hope ~ TREATMENT, df_TESS1)
model.Y <- lm(Ambition_norm ~ TREATMENT + hope + fear + anger, df_TESS1)
model.Y2 <- lm(Ambition_norm ~ TREATMENT + hope + fear + anger + party + AGE+ INCOME, df_TESS1)
results <- mediate(model.M, model.Y, treat='TREATMENT', mediator='hope', boot=TRUE, sims=1000, conf.level = .90, control.value = "Alexa", treat.value = "Clinton")
summary(results)
bt_effect <- c("ACME", "ADE", "Total Effect", 
               "Prop. Mediated")
bt_est <- c(results$d1, results$z1, results$tau.coef, results$n1)
#bt_p <- format.pval(c(fitMED$d1.p, fitMED$z1.p, fitMED$tau.p, fitMED$n1.p))
bt_p <- c(results$d1.p, results$z1.p, results$tau.p, results$n1.p)
bt_stars <- c(stars.pval(results$d1.p), stars.pval(results$z1.p),
              stars.pval(results$tau.p), stars.pval(results$n1.p))
bt_DF_Clinton <- data.frame(row.names = bt_effect, format(bt_est, digits = 2), 
                            format(bt_p, nsmall = 3), bt_stars)
colnames(bt_DF_Clinton) <- c("Coefficient", "p-value", "")

# Trump
df_TESS %>%
  filter(!is.na(hope), !is.na(Ambition_norm)) -> df_TESS1
model.0 <- lm(Ambition_norm ~ TREATMENT, df_TESS1)
model.M <- lm(hope ~ TREATMENT, df_TESS1)
model.Y <- lm(Ambition_norm ~ TREATMENT + hope + fear + anger, df_TESS1)
results <- mediate(model.M, model.Y, treat='TREATMENT', mediator='hope', boot=TRUE, sims=1000, conf.level = .90, control.value = "Alexa", treat.value = "Trump")
summary(results)
bt_effect <- c("ACME", "ADE", "Total Effect", 
               "Prop. Mediated")
bt_est <- c(results$d1, results$z1, results$tau.coef, results$n1)
bt_p <- c(results$d1.p, results$z1.p, results$tau.p, results$n1.p)
bt_stars <- c(stars.pval(results$d1.p), stars.pval(results$z1.p),
              stars.pval(results$tau.p), stars.pval(results$n1.p))
bt_DF_Trump <- data.frame(row.names = bt_effect, format(bt_est, digits = 2), 
                          format(bt_p, nsmall = 3), bt_stars)
colnames(bt_DF_Trump) <- c("Coefficient", "p-value", "")

# MeToo
results <- mediate(model.M, model.Y, treat='TREATMENT', mediator='hope', boot=TRUE, sims=1000, conf.level = .90, control.value = "Alexa", treat.value = "MeToo")
bt_effect <- c("ACME", "ADE", "Total Effect", 
               "Prop. Mediated")
bt_est <- c(results$d1, results$z1, results$tau.coef, results$n1)
bt_p <- c(results$d1.p, results$z1.p, results$tau.p, results$n1.p)
bt_stars <- c(stars.pval(results$d1.p), stars.pval(results$z1.p),
              stars.pval(results$tau.p), stars.pval(results$n1.p))
bt_DF_MeToo <- data.frame(row.names = bt_effect, format(bt_est, digits = 2), 
                          format(bt_p, nsmall = 3), bt_stars)
colnames(bt_DF_MeToo) <- c("Coefficient", "p-value", "")

# Marches
results <- mediate(model.M, model.Y, treat='TREATMENT', mediator='hope', boot=TRUE, sims=1000, conf.level = .90, control.value = "Alexa", treat.value = "March")
bt_effect <- c("ACME", "ADE", "Total Effect", 
               "Prop. Mediated")
bt_est <- c(results$d1, results$z1, results$tau.coef, results$n1)
bt_p <- c(results$d1.p, results$z1.p, results$tau.p, results$n1.p)
bt_stars <- c(stars.pval(results$d1.p), stars.pval(results$z1.p),
              stars.pval(results$tau.p), stars.pval(results$n1.p))
bt_DF_Marches <- data.frame(row.names = bt_effect, format(bt_est, digits = 2), 
                            format(bt_p, nsmall = 3), bt_stars)
colnames(bt_DF_Marches) <- c("Coefficient", "p-value", "")

print("TESS Results:")
bt_DF_Clinton
bt_DF_Trump
bt_DF_MeToo
bt_DF_Marches


##################### ANGER #####################
# Clinton:
model.0 <- lm(Ambition ~ TREATMENT, df_TESS)
model.M <- lm(anger ~ TREATMENT, df_TESS)
model.Y <- lm(Ambition_norm ~ TREATMENT + hope + fear + anger, df_TESS)
results <- mediate(model.M, model.Y, treat='TREATMENT', mediator='anger', boot=TRUE, sims=1000, conf.level = .90, control.value = "Alexa", treat.value = "Clinton")
bt_effect <- c("ACME", "ADE", "Total Effect", 
               "Prop. Mediated")
bt_est <- c(results$d1, results$z1, results$tau.coef, results$n1)
bt_p <- c(results$d1.p, results$z1.p, results$tau.p, results$n1.p)
bt_stars <- c(stars.pval(results$d1.p), stars.pval(results$z1.p),
              stars.pval(results$tau.p), stars.pval(results$n1.p))
bt_DF_Clinton <- data.frame(row.names = bt_effect, format(bt_est, digits = 2), 
                            format(bt_p, nsmall = 3), bt_stars)
colnames(bt_DF_Clinton) <- c("Coefficient", "p-value", "")

# Trump
results <- mediate(model.M, model.Y, treat='TREATMENT', mediator='anger', boot=TRUE, sims=1000, conf.level = .90, control.value = "Alexa", treat.value = "Trump")
bt_effect <- c("ACME", "ADE", "Total Effect", 
               "Prop. Mediated")
bt_est <- c(results$d1, results$z1, results$tau.coef, results$n1)
bt_p <- c(results$d1.p, results$z1.p, results$tau.p, results$n1.p)
bt_stars <- c(stars.pval(results$d1.p), stars.pval(results$z1.p),
              stars.pval(results$tau.p), stars.pval(results$n1.p))
bt_DF_Trump <- data.frame(row.names = bt_effect, format(bt_est, digits = 2), 
                          format(bt_p, nsmall = 3), bt_stars)
colnames(bt_DF_Trump) <- c("Coefficient", "p-value", "")

# MeToo
results <- mediate(model.M, model.Y, treat='TREATMENT', mediator='anger', boot=TRUE, sims=1000, conf.level = .90, control.value = "Alexa", treat.value = "MeToo")
bt_effect <- c("ACME", "ADE", "Total Effect", 
               "Prop. Mediated")
bt_est <- c(results$d1, results$z1, results$tau.coef, results$n1)
bt_p <- c(results$d1.p, results$z1.p, results$tau.p, results$n1.p)
bt_stars <- c(stars.pval(results$d1.p), stars.pval(results$z1.p),
              stars.pval(results$tau.p), stars.pval(results$n1.p))
bt_DF_MeToo <- data.frame(row.names = bt_effect, format(bt_est, digits = 2), 
                          format(bt_p, nsmall = 3), bt_stars)
colnames(bt_DF_MeToo) <- c("Coefficient", "p-value", "")

# Marches
results <- mediate(model.M, model.Y, treat='TREATMENT', mediator='anger', boot=TRUE, sims=1000, conf.level = .90, control.value = "Alexa", treat.value = "March")
bt_effect <- c("ACME", "ADE", "Total Effect", 
               "Prop. Mediated")
bt_est <- c(results$d1, results$z1, results$tau.coef, results$n1)
bt_p <- c(results$d1.p, results$z1.p, results$tau.p, results$n1.p)
bt_stars <- c(stars.pval(results$d1.p), stars.pval(results$z1.p),
              stars.pval(results$tau.p), stars.pval(results$n1.p))
bt_DF_Marches <- data.frame(row.names = bt_effect, format(bt_est, digits = 2), 
                            format(bt_p, nsmall = 3), bt_stars)
colnames(bt_DF_Marches) <- c("Coefficient", "p-value", "")

print("TESS Results:")
bt_DF_Clinton
bt_DF_Trump
bt_DF_MeToo
bt_DF_Marches

##################### FEAR #####################

# Clinton:
model.0 <- lm(Ambition_norm ~ TREATMENT, df_TESS)
model.M <- lm(fear ~ TREATMENT, df_TESS)
model.Y <- lm(Ambition_norm ~ TREATMENT + hope + fear + anger, df_TESS)
results <- mediate(model.M, model.Y, treat='TREATMENT', mediator='fear', boot=TRUE, sims=1000, conf.level = .90, control.value = "Alexa", treat.value = "Clinton")
bt_effect <- c("ACME", "ADE", "Total Effect", 
               "Prop. Mediated")
bt_est <- c(results$d1, results$z1, results$tau.coef, results$n1)
bt_p <- c(results$d1.p, results$z1.p, results$tau.p, results$n1.p)
bt_stars <- c(stars.pval(results$d1.p), stars.pval(results$z1.p),
              stars.pval(results$tau.p), stars.pval(results$n1.p))
bt_DF_Clinton <- data.frame(row.names = bt_effect, format(bt_est, digits = 2), 
                            format(bt_p, nsmall = 3), bt_stars)
colnames(bt_DF_Clinton) <- c("Coefficient", "p-value", "")

# Trump
results <- mediate(model.M, model.Y, treat='TREATMENT', mediator='fear', boot=TRUE, sims=1000, conf.level = .90, control.value = "Alexa", treat.value = "Trump")
bt_effect <- c("ACME", "ADE", "Total Effect", 
               "Prop. Mediated")
bt_est <- c(results$d1, results$z1, results$tau.coef, results$n1)
bt_p <- c(results$d1.p, results$z1.p, results$tau.p, results$n1.p)
bt_stars <- c(stars.pval(results$d1.p), stars.pval(results$z1.p),
              stars.pval(results$tau.p), stars.pval(results$n1.p))
bt_DF_Trump <- data.frame(row.names = bt_effect, format(bt_est, digits = 2), 
                          format(bt_p, nsmall = 3), bt_stars)
colnames(bt_DF_Trump) <- c("Coefficient", "p-value", "")

# MeToo
results <- mediate(model.M, model.Y, treat='TREATMENT', mediator='fear', boot=TRUE, sims=1000, conf.level = .90, control.value = "Alexa", treat.value = "MeToo")
bt_effect <- c("ACME", "ADE", "Total Effect", 
               "Prop. Mediated")
bt_est <- c(results$d1, results$z1, results$tau.coef, results$n1)
bt_p <- c(results$d1.p, results$z1.p, results$tau.p, results$n1.p)
bt_stars <- c(stars.pval(results$d1.p), stars.pval(results$z1.p),
              stars.pval(results$tau.p), stars.pval(results$n1.p))
bt_DF_MeToo <- data.frame(row.names = bt_effect, format(bt_est, digits = 2), 
                          format(bt_p, nsmall = 3), bt_stars)
colnames(bt_DF_MeToo) <- c("Coefficient", "p-value", "")

# Marches
results <- mediate(model.M, model.Y, treat='TREATMENT', mediator='fear', boot=TRUE, sims=1000, conf.level = .90, control.value = "Alexa", treat.value = "March")
bt_effect <- c("ACME", "ADE", "Total Effect", 
               "Prop. Mediated")
bt_est <- c(results$d1, results$z1, results$tau.coef, results$n1)
bt_p <- c(results$d1.p, results$z1.p, results$tau.p, results$n1.p)
bt_stars <- c(stars.pval(results$d1.p), stars.pval(results$z1.p),
              stars.pval(results$tau.p), stars.pval(results$n1.p))
bt_DF_Marches <- data.frame(row.names = bt_effect, format(bt_est, digits = 2), 
                            format(bt_p, nsmall = 3), bt_stars)
colnames(bt_DF_Marches) <- c("Coefficient", "p-value", "")

print("TESS Results:")
bt_DF_Clinton
bt_DF_Trump
bt_DF_MeToo
bt_DF_Marches

# ---------------------------------------------------------------------------------------------------------------
# ---------------------------------------------------------------------------------------------------------------

# Enthusiasm Averages on pages 21:
summarySE(df_TESS, measurevar="hope", groupvars=c("TREATMENT_sort"), conf.interval = .90)
t.test(x = df_TESS$hope[df_TESS$TREATMENT_sort == "Trump"], y = df_TESS$hope[df_TESS$TREATMENT_sort == "Alexa"], alternative = "less")
t.test(x = df_TESS$hope[df_TESS$TREATMENT_sort == "Clinton"], y = df_TESS$hope[df_TESS$TREATMENT_sort == "Alexa"])
t.test(x = df_TESS$hope[df_TESS$TREATMENT_sort == "March"], y = df_TESS$hope[df_TESS$TREATMENT_sort == "Alexa"])
t.test(x = df_TESS$hope[df_TESS$TREATMENT_sort == "MeToo"], y = df_TESS$hope[df_TESS$TREATMENT_sort == "Alexa"])

# ---------------------------------------------------------------------------------------------------------------

# Figure 4: Mean Enthusiasm by Experimental Condition
df_TESS1 <- subset(df_TESS,df_TESS$hope != 'NA') # Create a data frame for analysis
tgc <- summarySE(df_TESS1, measurevar="hope", groupvars=c("TREATMENT"), conf.interval = .90) 
tgc$low = tgc$hope-tgc$ci
tgc$high = tgc$hope+tgc$ci
tgc$sort <- c(5, 2, 3, 4, 1)
tgc %>%
  arrange(sort) -> tgc
# Generate Black and White Bar Plot
ggplot(data=tgc, aes(x=reorder(TREATMENT, sort), y=hope)) +
  geom_bar(position = "dodge", stat="identity") +
  geom_errorbar(position = "dodge",aes(ymin=low, ymax=high)) + 
  steph_theme_gray + 
  labs(title="Mean Levels of Enthusiasm", 
       x="Treatment Condition", y = "Average level of enthusiasm", caption = "Fall 2019 NORC Survey",subtitle = "By Experimental Condition", fill = "Gender:") + 
  theme(legend.position = "none")

# ---------------------------------------------------------------------------------------------------------------

# Table A4 referenced on page 21
Table_A4_A <- lm(hope ~ TREATMENT, data = df_TESS)
Table_A4_B <- lm(hope ~ TREATMENT*party, data = df_TESS)
export_summs(Table_A4_A,Table_A4_B, note= "+++ p < 0.001; ++ p < 0.01; + p < 0.05. one-tailed", stars = c(`+++` = 0.002, `++` = 0.02, `+` = 0.1)) #Set stars to double 0.05, 0.01, and 0.001 to show significance for one-tailed test rather than 2, note in the notes section.
export_summs(Table_A4_A,Table_A4_B, note= "*** p < 0.001; ** p < 0.01; * p < 0.05. two-tailed") # for two-tailed.

# ---------------------------------------------------------------------------------------------------------------

# Table A4b
df_TESS$hopeful_alone <- scale_norm(df_TESS$M1_4)
df_TESS$enthusiastic_alone <- scale_norm(df_TESS$M1_1) 
Table_4B_A <- lm(hopeful_alone~TREATMENT_sort, data = df_TESS)
Table_4B_B <- lm(enthusiastic_alone~TREATMENT_sort, data = df_TESS)
huxtable::huxreg(Table_4B_A,Table_4B_B,error_format = "({std.error})")
export_summs(Table_4B_A,Table_4B_B, note= "+++ p < 0.001; ++ p < 0.01; + p < 0.05. one-tailed", stars = c(`+++` = 0.002, `++` = 0.02, `+` = 0.1)) #Set stars to double 0.05, 0.01, and 0.001 to show significance for one-tailed test rather than 2, note in the notes section.
export_summs(Table_4B_A,Table_4B_B, note= "*** p < 0.001; ** p < 0.01; * p < 0.05. two-tailed") # for two-tailed
# ---------------------------------------------------------------------------------------------------------------

# Figure 5: Change in Enthusiasm Moving from the Control Group to the Clinton, Marches, #MeToo, and Trump Conditions, by Party
df_TESS1 <- subset(df_TESS, df_TESS$party != "NA") # Generate Data Frame for Analysis
hope1 <- lm(hope ~ TREATMENT*party, data = df_TESS1) # Generate linear model for HOPE/Enthusiasm
margins_summary(hope1, at = list(party = range(df_TESS1$party)))[3:10,] #Calculate Marginal Effects of Hope at 90% Confidence Interval
# Generate a dataframe with using the results of the marginal effects. We don't need to subdivide by gender because TESS
# was a sample of educated women.
data <- data.frame("model" = c("Republican Women","Democratic Women","Republican Women","Democratic Women","Republican Women","Democratic Women", "Republican Women","Democratic Women"), 
                   "estimate" = c(0.1877, 0.5071, 0.0724, 0.5256, 0.1049, 0.3180, -0.0844, -0.0724), 
                   "conf.low" = c(0.1135, 0.4514, 0.0001, 0.4697, 0.0325, 0.2623, -0.1613, -0.1292),
                   "conf.high" = c(0.2620, 0.5628, 0.1448, 0.5816, 0.1773, 0.3737, -0.0075, -0.0155), 
                   "sort" = c(1,2,3,4,5,6,7,8),
                   "term"=c("Clinton", "Clinton", "Marches","Marches", "#MeToo", "#MeToo", "Trump","Trump"))
# Plot a dot and whisker graph in Black and White for Enthusiasm.
dwplot(data) +
  theme_bw() +
  steph_theme_gray +
  theme(legend.justification=c(-1.0, .993), legend.position=c(-0.25, .99),
        legend.title = element_blank(), legend.background =
          element_rect(color="gray90")) + 
  xlab("") +
  scale_colour_grey(start = .3, end = .7,
                    name = "model") +
  labs(caption = "Fall 2019 NORC Survey") +
  geom_vline(xintercept = 0, colour = "grey60", linetype = 2) +
  ggtitle("Treatment Effect on Enthusiasm")

# ---------------------------------------------------------------------------------------------------------------

# Table A5 referenced on page 23
Table_A5 <- lm(anger ~ TREATMENT+party+AGE+INCOME, data = df_TESS)
export_summs(Table_A5, note= "+++ p < 0.001; ++ p < 0.01; + p < 0.05. one-tailed", stars = c(`***` = 0.002, `**` = 0.02, `*` = 0.1)) #Set stars to double 0.05, 0.01, and 0.001 to show significance for one-tailed test rather than 2, note in the notes section.
export_summs(Table_A5,error_format = "({std.error})") # two-tailed

# ---------------------------------------------------------------------------------------------------------------

# Table A6 as referenced on page 23
Table_A6 <- lm(hope ~ TREATMENT+party+AGE+INCOME, data = df_TESS)
export_summs(Table_A6,error_format = "({std.error})")
export_summs(Table_A6, note= "*** p < 0.001; ** p < 0.01; * p < 0.05. one-tailed", stars = c(`***` = 0.002, `**` = 0.02, `*` = 0.1)) #Set stars to double 0.05, 0.01, and 0.001 to show significance for one-tailed test rather than 2, note in the notes section.


# ---------------------------------------------------------------------------------------------------------------

# Table A7:
df_TESS %>%
  filter(!is.na(Ambition_norm)) -> df_TESS
# Do you think that what happens generally to women in this country will have something to do with what happens in your life?
# Will it affect you?
# Create the dummies for consistency:
df_TESS$CLINTON <- car::recode(df_TESS$TREATMENT, "'Alexa'= 0; 'Trump'=0; 'MeToo'=0; 'March'=0; 'Clinton'= 1")
df_TESS$TRUMP <- car::recode(df_TESS$TREATMENT, "'Alexa'= 0; 'Trump'=1; 'MeToo'=0; 'March'=0; 'Clinton'= 0")
df_TESS$METOO <- car::recode(df_TESS$TREATMENT, "'Alexa'= 0; 'Trump'=0; 'MeToo'=1; 'March'=0; 'Clinton'= 0")
df_TESS$MARCH <- car::recode(df_TESS$TREATMENT, "'Alexa'= 0; 'Trump'=0; 'MeToo'=0; 'March'=1; 'Clinton'= 0")
mix_mod <- '
Ambition_norm ~ b1 * Linked_Fate2B + c1 * TRUMP + c2 * CLINTON + c3 * MARCH + c4 * METOO
Linked_Fate2B ~ a1 * TRUMP + a99 * CLINTON + a98 * MARCH + a97 * METOO
indirTRUMPLinked := a1 * b1
indirCLINTONLinked := a99 * b1
indirMARCHLinked := a98 * b1
indirMETOOLinked := a97 * b1
# total effect
total := c1 + (a1 * b1)
'
fit <- lavaan:::sem(mix_mod, data=df_TESS, std.lv=TRUE)
summary(fit)

# ---------------------------------------------------------------------------------------------------------------

# Figure A1: Levels of Political Ambition by Treatment Condition Among TESS (NORC)
df_TESS1 <- subset(df_TESS,df_TESS$Ambition_norm != 'NA') # Create a data frame for analysis
tgc <- summarySE(df_TESS1, measurevar="Ambition_norm", groupvars=c("TREATMENT"), conf.interval = .90) 
tgc$low = tgc$Ambition_norm-tgc$ci
tgc$high = tgc$Ambition_norm+tgc$ci
tgc$sort <- c(5, 2, 3, 4, 1)
tgc %>%
  arrange(sort) -> tgc
ggplot(data=tgc, aes(x=reorder(TREATMENT, sort), y=Ambition_norm)) +
  geom_bar(position = "dodge", stat="identity") +
  geom_errorbar(position = "dodge",aes(ymin=low, ymax=high)) + 
  steph_theme_gray + 
  labs(title="Mean Levels of Political Ambition", 
       x="Treatment Condition", y = "Average Level of Political Ambition", caption = "Fall 2019 NORC Survey",subtitle = "By Experimental Condition", fill = "Gender:") + 
  theme(legend.position = "none")

# -----------------------------------------------------------------------------------------------------------

# Load in the CCES 2020 Data:
df <- readRDS("CCES2020_Data.rds")


 # Table A8. Modeling non-compliance by treatment, gender, and party
df %>%
  filter(!is.na(female))-> df_sub
mod1 <- lm(noncomplier ~ TREATMENT * female, data = df_sub, weights = teamweight)
df %>%
  filter(!is.na(dem)) -> df_sub
mod2 <- lm(noncomplier ~ TREATMENT * dem, data = df_sub, weights = teamweight)
# Both?
df %>%
  filter(!is.na(dem),
         !is.na(female)) -> df_sub
mod3 <- lm(noncomplier ~ TREATMENT * dem * female, data = df_sub, weights = teamweight)
export_summs(mod1,mod2,mod3, note= "+++ p < 0.001; ++ p < 0.01; + p < 0.05. one-tailed", stars = c(`+++` = 0.002, `++` = 0.02, `+` = 0.1)) #Set stars to double 0.05, 0.01, and 0.001 to show significance for one-tailed test rather than 2, note in the notes section.
export_summs(mod1,mod2,mod3, note= "*** p < 0.001; ** p < 0.01; * p < 0.05. two-tailed", stars = c(`***` = 0.001, `**` = 0.01, `*` = 0.05)) # two-tailed

# -----------------------------------------------------------------------------------------------------------

# Table A9. Average marginal effects (AME) on noncompliance by treatment and party
df$dem <- as.numeric(as.character(df$dem))
df %>%
  filter(!is.na(dem)) -> df_sub
mod2 <- lm(noncomplier ~ TREATMENT * dem, data = df_sub, weights = teamweight)
margins_summary(mod2, at = list(dem = range(df_sub$dem)))[4:8,]

# -----------------------------------------------------------------------------------------------------------

# Table A10. Average marginal effects (AME) on noncompliance by treatment and gender
df %>%
  filter(!is.na(female)) -> df_sub
mod1 <- lm(noncomplier ~ TREATMENT * female, data = df_sub, weights = teamweight)
margins_summary(mod1, at = list(female = range(df_sub$female)))[4:8,]

# -----------------------------------------------------------------------------------------------------------

#Table A11. Average marginal effects (AME) on noncompliance by treatment, gender, and party

df %>%
  filter(!is.na(female), !is.na(dem)) -> df_sub
mod1 <- lm(noncomplier ~ TREATMENT * female * dem, data = df_sub, weights = teamweight)
margins_summary(mod1, at = list(female = range(df_sub$female), dem = range(df_sub$dem)))[9:20,]

# -----------------------------------------------------------------------------------------------------------

# Table A12. CCES anger induction on ambition (Two-Stage Least Squares Regression)
df %>%
  mutate(
    UCR121_treatment = as.factor(df$TREATMENT),
    UCR121_treatmentc = as.factor(df$TREATMENT),
    Relaxed = if_else(UCR121_treatment == 'Relaxed', 1, 0),
    Treatment_of_Minorities = if_else(UCR121_treatment == 'Treatment of Minorities', 1, 0),
    Treatment_of_Women = if_else(UCR121_treatment == 'Treatment of Women', 1, 0),
    Trump = if_else(UCR121_treatment == 'Trump', 1, 0),
    Relaxedc = if_else(noncomplier == 1, Relaxed, 0),
    Treatment_of_Minoritiesc = if_else(noncomplier == 1, Treatment_of_Minorities, 0),
    Treatment_of_Womenc = if_else(noncomplier == 1, Treatment_of_Women, 0),
    Trumpc = if_else(noncomplier == 1, Trump, 0)
  ) -> df

mod1 <- iv_robust(ambition ~ Treatment_of_Minoritiesc + Treatment_of_Womenc + Trumpc | Treatment_of_Minorities + Treatment_of_Women + Trump, data = df, se_type = "stata", weights = teamweight)
export_summs(mod1, note= "+++ p < 0.001; ++ p < 0.01; + p < 0.05. one-tailed", stars = c(`+++` = 0.002, `++` = 0.02, `+` = 0.1)) #Set stars to double 0.05, 0.01, and 0.001 to show significance for one-tailed test rather than 2, note in the notes section.
export_summs(mod1, note= "*** p < 0.001; ** p < 0.01; * p < 0.05. two-tailed", stars = c(`***` = 0.001, `**` = 0.01, `*` = 0.05)) # two-tailed
